Sea $A \in \real^{n \times n}$, se desea hallar una descomposici\'on de 
la matr\'iz $A$ en tres matrices, si existen, 
$L,U,P \in \times^{n \times n}$ tales que $L$ es triangular inferior, 
$U$ es triangular superior y $P$ es una matriz de permutacion. 
Esta descomposici\'on es tal que $PA = LU$. 

\paragraph{}
\textbf{Eliminaci\'on Gaussiana}
\par
Un m\'etodo para resolver sistemas de ecuaciones lineales es el 
m\'etodo de eliminaci\'on de Gauss. 
La idea es manipular las ecucaciones por medio de operaciones 
elementales para transformar el sistema original en un sistema 
equivalente que sea m\'as sencillo de resolver (sistema triangular). 
En este caso se operar\'a con funciones de $\real^n \mapsto \real^n$, 
es decir, tales que el sistema matricial resulte con una matr\'iz 
$A$ cuadrada.

\paragraph{}
Las operaciones elementales en la eliminaci\'on de Gauss son tres:
\begin{enumerate}
    \item Multiplicaci\'on de una ecuación por una constante no nula.
    \item Sustracción del m\'ultiplo de otra ecuaci\'on a la ecuaci\'on 
			en cuesti\'on.
    \item Intercambio de ecuaciones.
\end{enumerate}

Estas mismas operaciones elementales pueden traducirse a la forma 
matricial, realizando en este caso operaciones sobre las filas de las 
matrices.
El sistema inicial es $A^{(1)}x$ $=$ $b^{(1)}$. Los supra\'indices se 
utilizan para indicar el estado del sistema.

%~ \paragraph{}
%~ \textbf{Primer paso de eliminaci\'on}
%~ Si $\displaystyle a_{11}^{(1)} \neq 0$, se puede eliminar la inc\'ognita $x_1$ de las dem\'as ecuaciones. El paso t\'ipico es restar de la $i$-\'esima ecuación ($i = 2, 3, \dots, n$) la primera multiplicada por
%~ \begin{equation}
%~ \displaystyle m_{i1} =  \frac{a_{i1}^{(1)}}{a_{11}^{(1)}}  \quad \forall \; i \quad i = 2, 3, \dots, n
%~ \end{equation}
%~ 
%~ Luego de estas restas, la $i$-\'esima ecuaci\'on tendr\'a nuevos coeficiente $a_{ij}^{(2)}$ y $b_{i}^{(2)}$ cuyos valores son:
%~ \begin{align}
%~ a_{i1}^{(2)} &= 0 \quad \forall \; i \quad i = 2, 3, \dots, n \\
%~ a_{ij}^{(2)} &= a_{ij}^{(1)} - m_{i1} . a_{1j}^{(1)} \quad \forall \; i, j \quad i,j = 2, 3, \dots, n\\
%~ b_{i}^{(2)} &= b_{i}^{(1)} - m_{i1} . b_{1}^{(1)} \quad \forall \; i \quad i = 2, 3, \dots, n
%~ \end{align}
%~ 
%~ Se obtiene entonces el nuevo sistema $A^{(2)}x$ $=$ $b^{(2)}$
%~ 
%~ \begin{center}
%~ $\begin{array}{rcl}
     %~ a_{11}^{(1)}.x_1 + a_{12}^{(1)}.x_2 +  \ldots  + a_{1n}^{(1)}.x_n & = & b_1^{(1)} \\
                        %~ a_{22}^{(2)}.x_2 +  \ldots  + a_{2n}^{(2)}.x_n & = & b_2^{(2)} \\
                            %~ \vdots  & \vdots & \vdots \\
                        %~ a_{n2}^{(2)}.x_2 +  \ldots + a_{nn}^{(2)}.x_n & = & b_n^{(2)}
%~ \end{array}$
%~ \end{center}

%~ \paragraph{}
%~ \textbf{Segundo paso de eliminación}
%~ En este paso el objetivo es eliminar la inc\'ognita $x_2$ de la tercera a la \'ultima ecuación. Si $a_{22}$ $\neq$ $0$, se resta de la $i$-\'esima ecuación ($i = 3, \dots, n$) la segunda multiplicada por
%~ \begin{equation}
%~ m_{i2} =  a_{i2}^{(2)} / a_{22}^{(2)}  \quad \forall \; i \quad i = 3, \dots, n
%~ \end{equation}
%~ 
%~ Los nuevos coeficientes $a_{ij}^{(3)}$ y $n_{i}^{(3)}$ de la $i$-ésima ecuación serán
%~ \begin{align}
%~ a_{i2}^{(3)} &= 0 \quad \forall \; i \quad i = 3, \dots, n \\
%~ a_{ij}^{(3)} &= a_{ij}^{(2)} - m_{i2} . a_{2j}^{(2)} \quad \forall \; i, j \quad i,j = 3, \dots, n\\
%~ b_{i}^{(3)} &= b_{i}^{(2)} - m_{i2} . b_{2}^{(2)} \quad \forall \; i \quad i = 3, \dots, n
%~ \end{align}
%~ 
%~ Se obtiene entonces el nuevo sistema $A^{(3)}x$ $=$ $b^{(3)}$

%~ \paragraph{}
%~ \textbf{Siguientes pasos de eliminaci\'on}
%~ Repitiendo lo anterior, se van obteniendo los sistemas \\$A^{(i)}x$ = $b^{(i)}$ y cuando $i = n$ (es decir despu\'es de $n-1$ pasos de eliminaci\'on) se obtiene un sistema triangular superior.

%~ \begin{center}
%~ $
%~ \displaystyle
%~ \begin{array}{rcl}
%~ a_{11}^{(1)}.x_1 + a_{12}^{(1)}.x_2 + a_{13}^{(1)}.x_3 +  \ldots  + a_{1n}^{(1)}.x_n & = & b_1^{(1)} \\
                    %~ a_{22}^{(2)}.x_2 + a_{23}^{(2)}.x_3 +  \ldots  + a_{2n}^{(2)}.x_n & = & b_2^{(2)} \\
                                      %~ + a_{33}^{(3)}.x_3 +  \ldots  + a_{3n}^{(3)}.x_n & = & b_3^{(3)} \\
                                            %~ \ddots & \vdots & \vdots \\
                                                                        %~ a_{nn}^{(n)}.x_n & = & b_n^{(n)}
%~ \end{array} $
%~ \end{center}

El proceso de eliminaci\'on termina sin problemas siempre que ninguno 
de los pivotes $a_{11}^{(1)}$, $a_{22}^{(2)}$, \dots , $a_{nn}^{(n)}$ 
sea cero. 
Este sistema triangular superior $Ux$ = $b$ (donde $U$ = $A^{(n)}$, 
$b$ = $b^{(n)}$ ) equivalente al sistema original. 
Sin embargo, este nuevo sistema puede resolverse muy fácilmente por 
medio de la técnica de sustitución hacia atrás (backward substitution):

\begin{align}
	x_n &= \displaystyle \frac{b_n}{a_{nn}} \\
	x_i &= \displaystyle \frac{b_i - \displaystyle 
			\sum^{n}_{\substack{k=i+1}} a_{ij}x_j}{a_{ii}}, \quad 
			\forall \; i \quad i = n-1, n-2 \dots, 1
\end{align}

\paragraph{}
La complejidad de resolver el sistema $\displaystyle Ax = b$ es de 
$O(n^{3})$. 
Al obtener la factorizacion $PLU$ de $A$ es posible resolver el sistema 
equivalente en $O(n^{2})$. 
%~ \par
Esto se realiza tomando $\displaystyle y = Ux$ y resolviendo los 
sistemas $\displaystyle Ly = Pb$ y $\displaystyle Ux=y$ ambos en 
$O(n^{2})$. 
Este hecho se debe a que $U$ y $L$ son matrices triangulares.
 %~ (el primer sistema se resuelve con \textsl{forward substitution} y el segundo con \textsl{backward substitution}).
\par
Si se deben resolver m\'ultiples instancias de la forma $Ax=b_{i}$, 
el costo ser\'a de $O(n^{3})$ para la primera (para obtener la 
descomposici\'on), y $O(n^{2})$ para las siguientes instancias. 

\paragraph{}
Si es posible efectuar la eliminaci\'on gaussiana para el sistema 
$Ax = b$ sin intercambio de filas, entonces $A$ es factorizable como el 
producto de $L$ y $U$. 
Si se requiere intercambio de fila, puede hallarse $P$ (matr\'iz de 
permutaci\'on) tal que $PA = LU$. 
Para garantizar que esta factorizaci\'on sea \'unica por lo general se 
pide que haya unos en la diagonal de L (Doolitle) o en la de U (Crout).

\paragraph{}
Durante el algoritmo de eliminacion gaussiana, pueden ir 
almacen\'andose los coeficientes usados para obtener ceros en una 
columna en el mismo lugar donde deber\'ia haber un nuevo cero, 
logrando de esta forma un uso eficiente de la memoria (en una sola 
matr\'iz se tiene la informaci\'on de $L$ y de $U$).

\begin{theorem}
	Una matriz A no singular tiene factorizaci\'on $LU$ si cumple 
	alguna de las siguientes condiciones equivalentes equivalencias:
	\begin{itemize}
		\setlength{\itemsep}{0pt}
		\item{No emerge un pivote igual a cero durante la 
				eliminaci\'on gaussiana.}
		\item{Todas las submatrices principales de $A$ son 
				no singulares.}
	\end{itemize}
\end{theorem}

\begin{remark}
	Para toda matriz $A$ no singular existe una matriz de permutaci\'on 
	$P$ tal que $PA$ tiene \\factorizaci\'on $LU$.
\end{remark}

La matr\'iz $U$ es de la forma 
$U=M_{(n-1)}M_{(n-2)}\dots M_{(2)}M_{(1)}A_{(1)}$,
donde $M_{k}$ es la matr\'iz utilizada en el k-\'esimo paso de la 
eliminaci\'on gaussiana. 
$M_{k}$ es inversible por ser producto de matrices elementales (las 
matrices utilizadas para realizar las tres operaciones elementales 
sobre matrices), tiene unos en la diagonal, todos ceros excepto en la 
columna $k$ que contiene, debajo de la diagonal, los multiplicadores 
(negados) usados para colocar ceros en cada fila de $A_{(k)}$.

La matr\'iz $L$ es de la forma 
${M_{(1)}}^{-1} {M_{(2)}}^{-1} \dots {M_{(n-2)}}^{-1} {M_{(n-1)}}^{-1}$,
que resulta ser triangular inferior con unos en su diagonal.
cuyos \'unicos elementos debajo de la diagonal son los multiplicadores 
usados para triangular $A$.

\paragraph{}
Este m\'etodo sufre de inestabilidad num\'erica cuando la matriz $A$
tiene un numero de condici\'on grande. Se puede utilizar pivoteo parcial
(intercambiar filas para dividir por el numero de mayor modulo en la 
primera columna de la submatr\'iz a factorizar), o pivoteo total 
(intercambiar filas y columnas para dividir por el numero de mayor 
múdulo ubicado en la submatríz a factorizar), para mejorar la 
estabilidad num\'erica de este m\'etodo (que no es m\'as que la 
eliminaci\'on gaussiana). 

%~ La estrategia de pivoteo parcial consiste en elegir en cada paso como pivote el elemento de la columna que resulte mayor en m\'odulo.
%~ 
%~ \begin{enumerate}    
    %~ \item Para cada columna $j$ 
    %~ \begin{enumerate}
        %~ \item Elegir como pivote el elemento de la columna que sea mayor en m\'odulo.\\
        %~ pivote[$j$] $=$ $p$ $=$ $\displaystyle \max_{i = j, \dots, n}{\displaystyle |a_{ij}|}$
        %~ \item Intercambiar fila $j$ con fila $p$ en la matr\'iz $A$.
        %~ \item Intercambiar fila $j$ con fila $p$ en $b$.
        %~ \item Realizar el paso de eliminaci\'on gaussiana normalmente ahora que el pivote deseado esta ubicado en $A[j][j]$.
    %~ \end{enumerate}    
%~ \end{enumerate}
